/**
This file is part of Gel2DE.

    Gel2DE is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Gel2DE is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with Gel2DE.  If not, see <http://www.gnu.org/licenses/>.

	Copyright (C) 2012, University of Bergen
*/


#include "SpearmanRankOrderCorrCoeff.h"
#include <itkTDistribution.h>

SpearmanRankOrderCorrCoeff::SpearmanRankOrderCorrCoeff(int aNelements)
{
	m_data = new SpearmanRankData(aNelements);	
}

SpearmanRankOrderCorrCoeff::SpearmanRankOrderCorrCoeff(std::vector<double>* aX, std::vector<double>* aY)
{
	m_data = new SpearmanRankData(aX, aY);
}
SpearmanRankOrderCorrCoeff::~SpearmanRankOrderCorrCoeff(void)
{
	delete m_data;
}

void SpearmanRankOrderCorrCoeff::Clear()
{
	m_data->Clear();
}

//Method for calculating r_s (correlation coefficient)
//and prob_rs (significance).
//returns a pair of <r_s, prob_rs>
std::pair<double, double>* SpearmanRankOrderCorrCoeff::calculateCorrCoeff(bool aCalcSignificance)
{
	m_data->RankData();
	//Ready to compute
	std::vector<SpearmanRankElement*>::iterator i = m_data->GetData()->begin();
	double D = 0;
	for(i; i != m_data->GetData()->end(); ++i)
	{
		double dRank = (*i)->GetRank(SpearmanRankElement::X) - (*i)->GetRank(SpearmanRankElement::Y);
		D += dRank * dRank;
	}
		
	double N = m_data->GetData()->size();
	double N3N = N*N*N - N;

	double sumf3f = 0;
	std::vector<int>* m_nbrXties = m_data->GetXties();
	std::vector<int>::iterator iterf = m_nbrXties->begin();
	for(iterf ; iterf != m_nbrXties->end(); iterf++)
	{
		double f3f = (*iterf)*(*iterf)*(*iterf) - (*iterf);
		sumf3f += f3f;			
	}
	double sumg3g = 0;
	std::vector<int>* m_nbrYties = m_data->GetYties();
	std::vector<int>::iterator iterg = m_nbrYties->begin();
	for(iterg ; iterg != m_nbrYties->end(); iterg++)
	{
		double g3g = (*iterg)*(*iterg)*(*iterg) - (*iterg);
		sumg3g += g3g;			
	}

	double top = D + (sumf3f + sumg3g)/12.;
	top *= 6./N3N;
	top = 1 - top;

	double bottomA = sqrt(1 - (sumf3f/N3N));
	double bottomB = sqrt(1 - (sumg3g/N3N));
	double bottom = bottomA * bottomB;
	double rs = 0;
	if(bottom != 0.)
	{
		rs = top/bottom;
	}

	double t;
	double probrs;
	if(fabs(rs)<1.0)
	{
		t = rs*sqrt(((double)N - 2.)/(1. - (rs*rs)));

		if (aCalcSignificance)
		{

			if(rs < 0)
				probrs = 2*itk::Statistics::TDistribution::CDF(t, N - 2.);
			else
				probrs = 2* (1. - itk::Statistics::TDistribution::CDF(t, N - 2.));
		}
	}
	else
	{
		probrs = 0;
	}
			
	return new std::pair<double, double>(rs, probrs);

}

void SpearmanRankOrderCorrCoeff::AddData(std::vector<double>* aX, std::vector<double>* aY)
{
	m_data->AddElements(aX, aY);
}

void SpearmanRankOrderCorrCoeff::AddData(double aX, double aY)
{
	m_data->AddElement(aX, aY);
}


double SpearmanRankOrderCorrCoeff::StdDevI()
{
	//todo: reduce from 2N to N.
	std::vector<SpearmanRankElement*>* data = m_data->GetData();
	std::vector<SpearmanRankElement*>::iterator iter = data->begin();
	double sum = 0;
	for(iter; iter != data->end(); iter++)
	{
		sum += (*iter)->GetValue(SpearmanRankElement::X);
	}
	int N = data->size();
	double avg = sum/N;
	if(avg == 0) return 0;

	double devSum = 0;
	for(iter = data->begin(); iter != data->end(); iter++)
	{
		double a = (*iter)->GetValue(SpearmanRankElement::X) - avg;
		devSum += a*a;
	}

	double stdDev = sqrt(devSum/(N - 1));
	//std::cout << " stdDev : " << stdDev << std::endl;
	return stdDev;//*avg;
}


